clear all;clc
tic

er0=8.854187817e-12;
mu0=4*pi*1e-7;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NodePt=dlmread('data/node.txt');
PlaneNd=dlmread('data/plane.txt');
SurfaceNd=dlmread('data/surface.txt');

PlaneSize=size(PlaneNd);
SurfaceSize=size(SurfaceNd);

for n=1:PlaneSize(1)
    for ni=2:4
        PlaneNd(n,ni)=PlaneNd(n,ni)+1;
    end
end

for n=1:SurfaceSize(1)
    for ni=2:3
        SurfaceNd(n,ni)=SurfaceNd(n,ni)+1;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure
% hold on
% for n=1:PlaneSize(1)
%     for np=1:3
%         nodex(np)=NodePt(PlaneNd(n,np+1),2);
%         nodey(np)=NodePt(PlaneNd(n,np+1),3);
%     end     
%     plot([nodex(1) nodex(2) nodex(3) nodex(1)],[nodey(1) nodey(2) nodey(3) nodey(1)])
% end
% 
% for n=1:SurfaceSize(1) 
%     for np=1:2
%         nodex(np)=NodePt(SurfaceNd(n,np+1),2);
%         nodey(np)=NodePt(SurfaceNd(n,np+1),3);
%     end
%     plot([nodex(1) nodex(2) nodex(1)],[nodey(1) nodey(2) nodey(1)],'r','linewidth',2)
% end
% axis equal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[GlobKe,GlobBe]=CalculatePlaneKeBe(NodePt,PlaneNd);
[GlobKs,GlobBs]=CalculateSurfaceKeBe(NodePt,SurfaceNd);
GlobKe=GlobKe+GlobKs;
GlobBe=GlobBe+GlobBs;
phi=GlobKe\GlobBe;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
tri0=delaunay(NodePt(:,2),NodePt(:,3));
trisurf(tri0,NodePt(:,2),NodePt(:,3),phi)
%shading interp
view([0 0 1])
alpha(0.6)
axis equal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
toc